Les CPK

La papier “Discovery of nitrate–CPK–NLP signalling in central nutrient–growth networks” met en avant la découverte des CPK, des Calcium (rôle de signalement cellulaire) sensor Protein Kinases (molécule dont le role est de modifier chimiquement d’autres molécules cibles ajout l’ajout d’un groupe phosphate). Elles jouent un role de master regulators dans la réponse primaire au nitrate. En particulier, elles vont phosphoryler les facteurs de transcription NLP qui orchestrent de nombreux processus de nutrition et de développement.

Les auteurs créent des mutants CPK10, CPK30, CPK32, et NLP7 pour le prouver.

Préparation du jeu de données

On loade la matrice de reads par échantillon et par gènes obtenue par htseq-count. On supprime les dernières lignes qui correspondent aux reads non identifiés comme des gènes ou comme ambigus.

data = read.table("quantif_all.txt", h = T)
rownames(data) = data$Gene
data = dplyr::select(data, -c("Gene", "WT.K1", "WT.K2", "WT.K3"))
genes = which(!(grepl("__", rownames(data))))
data = data[genes, ]
data = data[, sort(colnames(data))]


data_mutant = read.table("quantif_mutants.txt", h = F)
colnames(data_mutant) <- c("Gene", "icpk.N1", "icpk.N2", "icpk.N3")
rownames(data_mutant) = data_mutant$Gene
data_mutant = dplyr::select(data_mutant, -c("Gene"))
data_mutant = data_mutant[genes, ]

for (var in colnames(data_mutant)) {
    data[, var] = data_mutant[, var]
}

kable(head(data))

Normalisation

In the biological point of view, a gene must be expressed at some minimal level before it is likelyto be translated into a protein or to be biologically important. In addition, the pronounceddiscreteness of these counts interferes with some of the statistical approximations that areused later in the pipeline. These genes should be filtered out prior to further analysis.As a rule of thumb, genes are dropped if they can’t possibly be expressed in all the samplesfor any of the conditions. On met aussi bien à jour la nouvelle taille de librairie.

y <- DGEList(counts = data, group = c(1, 1, 1, 2, 2, 2))
y$samples

        group lib.size norm.factors
WT.N1       1   406584            1
WT.N2       1   400980            1
WT.N3       1   493383            1
icpk.N1     2   376432            1
icpk.N2     2   472267            1
icpk.N3     2   557570            1

keep <- filterByExpr(y)
y <- y[keep, , keep.lib.sizes = FALSE]

RNA-seq provides a measure of the relative abundance of each gene in each RNAsample, but does not provide any measure of the total RNA output on a per-cell basis. This commonly becomes important when a small number of genes are very highly expressed in one sample, but not in another. The highly expressed genes can consume a substantial proportionof the total library size, causing the remaining genes to be under-sampled in that sample.

The calcNormFactorsfunction normalizes for RNA composition by finding a set of scalingfactors for the library sizes that minimize the log-fold changes between the samples for mostgenes. The default method for computing these scale factors uses a trimmed mean of M-values (TMM) between each pair of samples [30]. We call the product of the original library size and the scaling factor the effective library size. The effective library size replaces theoriginal library size in all downsteam analyses. TMM is recommended for most RNA-Seq data where the majority (more than half) of thegenes are believed not differentially expressed between any pair of the sample

y <- calcNormFactors(y)
y$samples

        group lib.size norm.factors
WT.N1       1   373992    1.0217476
WT.N2       1   370605    0.9818326
WT.N3       1   456687    0.9680863
icpk.N1     2   347890    1.0147058
icpk.N2     2   439399    0.9706227
icpk.N3     2   514209    1.0454764

norm <- cpm(y, normalized.lib.sizes = TRUE)
not_norm <- cpm(y, normalized.lib.sizes = FALSE)
rd_genes = sample(rownames(norm), 500)
heatmap(not_norm[rd_genes, ], main = "Sans normalisation")

heatmap(norm[rd_genes, ], main = "Avec normalisation")

The normalization factors of all the libraries multiply to unity. A normalization factor belowone indicates that a small number of high count genes are monopolizing the sequencing,causing the counts for other genes to be lower than would be usual given the library size. Asa result, the library size will be scaled down, analogous to scaling the counts upwards in thatlibrary. Conversely, a factor above one scales up the library size, analogous to downscalingthe counts.

The GC-content of each gene does not change from sample to sample, so it can be expectedto have little effect on differential expression analyses to a first approximation. Recent pub-lications, however, have demonstrated that sample-specific effects for GC-content can bedetected [28, 13]. The EDASeq [28] and cqn [13] packages estimate correction factors thatadjust for sample-specific GC-content effects in a way that is compatible with edgeR.

Like GC-content, gene length does not change from sample to sample, so it can be expectedto have little effect on differential expression analyses. Nevertheless, sample-specific effectsfor gene length have been detected [13], although the evidence is not as strong as for GC-content.

Estimation de la variance et Expression différentielle

Biological Variance is the coefficient of variation with which the (unknown) true abundanceof the gene varies between replicate RNA samples.

The dispersion estimates can be viewed in a BCV plot:

Estimation compliquée à documenter

y <- estimateDisp(y)

Design matrix not provided. Switch to the classic mode.

et <- exactTest(y)
best = topTags(et, n = 2000)
best = best[best$table$FDR < 0.05, ]
DEgenes = rownames(best$table)
plotBCV(y)

plotMD(et)
abline(h = c(-1, 1), col = "blue")
abline(h = c(0), col = "red")

summary(decideTests(et))

        2-1
Down    434
NotSig 5736
Up      369

On a 434+369 gènes différentiellement exprimés avec un seuil pour la pvalue ajustée de 5%.

The functionplotMDS draws a multi-dimensional scaling plot of the RNA samples in whichdistances correspond to leading log-fold-changes between each pair of RNA samples. Theleading log-fold-change is the average (root-mean-square) of the largest absolute log-fold-changes between each pair of sample

plotMDS(y, top = 100, col = rep(1:2, each = 3))

logcpm <- cpm(y, log = TRUE)
heatmap(logcpm[DEgenes, ])

On peut maintenant chercher à voir les ontologies des gènes sélectionnés :

library(knitr)
mart = useMart(biomart = "plants_mart", host = "plants.ensembl.org", dataset = "athaliana_eg_gene")
results <- getBM(filters = "ensembl_transcript_id", attributes = c("ensembl_transcript_id", 
    "description", "external_gene_name"), values = DEgenes, mart = mart)
rownames(results) = results$ensembl_transcript_id
r = results[DEgenes, ]
r$pval = best$table$PValue
kable(r, caption = "Genes DE et leur GO")

library(stringi)
DEGenesNames = stri_remove_empty(r$external_gene_name, na_empty = FALSE)
save(file = "DEGenesEdgeRMutant.RData", DEGenesNames)

Sur ces 803 gènes DE, 30 contiennent le mot “kinase” dans leur description, 8 le mot calcium.

library(data.table)
kable(r[r$description %like% "kinase", ])

kable(r[r$description %like% "Calcium" | r$description %like% "calcium", ])

kable(r[r$external_gene_name %like% "NLP", ])